Synchronization of oscillators with long range power law interactions 
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We present analytical calculations and numerical simulations for the synchronization of oscillators 
interacting via a long range power law interaction on a one dimensional lattice. We have identified 
the critical value of the power law exponent a c across which a transition from a synchronized to an 
unsynchronized state takes place for a sufficiently strong but finite coupling strength in the large 
system limit. We find a c = 3/2. Frequency entrainment and phase ordering are discussed as a 
function of a > 1. The calculations are performed using an expansion about the aligned phase state 
(spin-wave approximation) and a coarse graining approach. We also generalize the spin-wave results 
to the rf-dimensional problem. 
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I. INTRODUCTION 

Synchronization between a collection of oscillating ob- 
jects is a common feature in a number of complex sys- 
tems, such as the pacemaker cells in the heart, neurons in 
the nervous system, an array of Joscphson junctions and 
rhythmic applause in a theater [1]. At the same time, 
from a purely theoretical point of view, the phenomenon 
is interesting as perhaps the simplest examples of a col- 
lective response of driven dynamical systems, showing 
many features reminiscent of equilibrium phase transi- 
tions. Recent technological advancements in the fabrica- 
tion of nanoelectromcchanical systems (NEMS) promise 
large arrays of interacting nonlinear oscillators [2] that 
will provide good testing grounds for the theoretical pre- 
dictions, as well as potential applications in sensitive de- 
tectors and precise frequency sources. 

Winfree [3] introduced a simple phase description of 
coupled oscillators, and he and Kuramoto [4] demon- 
strated synchronization for the case in which each os- 
cillator is coupled equally to all the other oscillators (all- 
to-all coupling) using a mean field theory (for a review, 
see ref . [5] ) : above a critical coupling strength depending 
on the distribution of the oscillator frequencies a finite 
fraction of the oscillators become entrained and oscillate 
at the same frequency, leading to a coherent signal in the 
summed response of the oscillators. On the other hand, 
for nearest neighbor coupling Daido [6] and Strogatz and 
Mirollo [7] showed that there is no macroscopic entrain- 
ment for a one dimensional chain. (By macroscopic we 
mean an O(N) value for N — > oo oscillators.) In higher 
dimensions Strogatz and Mirollo further showed that any 
macroscopic entrained cluster must have the form of a 
sponge, i.e. any compact macroscopic entrained region 
contains "holes" of unentraincd oscillators. The full be- 
havior of the nearest neighbor model as a function of the 
dimension of the lattice is not completely understood, al- 
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though based on approximate analytic arguments [6, 8] 
and numerics [9] it is conjectured that d = 2 is the lower 
critical dimension for which macroscopic entrainment oc- 
curs. 

We consider the phase model for a one-dimensional 
chain of oscillators with a strength of coupling between 
oscillators that falls off as a power law of their separation, 
and study how the synchronization behavior changes as 
a function of the power. This system was introduced by 
Rogers and Wille [10]. They investigated the system nu- 
merically, in particular finding the the critical value of the 
power law for synchronization of a macroscopic system 
as a function of the coupling strength. They suggested a 
critical value a c = 2 above which macroscopic synchro- 
nization would not occur for any finite coupling strength, 
no matter how large. Marodi ct al. [11] further investi- 
gated the system, focusing particularly on the question of 
complete entrainment in finite systems for small enough 
power laws. For an interaction power law a < 1, the 
sum of the interactions of a single oscillator with an in- 
finite lattice of oscillators with aligned phases diverges 
for fixed coupling strength. In their study of this range 
of a, Rogers and Wille [10] chose a system size depen- 
dent normalization of the coupling strength to remove 
this divergence. This choice of normalization allows the 
investigation of the crossover to the all-to-all model of N 
oscillators, where the coupling to each oscillator is scaled 
by N —1 so that the synchronization transition occurs at 
finite coupling constant. On the other hand, Marodi et 
al. [11] argued that in physical systems of interest, the 
interaction strength would not be expected to scale with 
the system size, and they investigated the model without 
scaling the coupling constant with system size A major 
focus of their study was how the size of the system needed 
for complete entrainment then depends on the interaction 
power law and the coupling strength. 

A power law interaction is of interest both exper- 
imentally and theoretically. This type of interaction 
should be relevant to some implementations of nanome- 
chanical arrays, since both electrostatic interactions be- 
tween charges or dipoles on the devices and elastic in- 
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teractions through the supporting substrate may lead to 
such long range interactions. Interactions falling off as 
a power law have also been used to model the complex 
long range connectivity of neurons [12]. Radicchi and 
Meyer-Ortmanns [13] have recently investigated the case 
of a single pacemaker oscillator with a different frequency 
coupled through a power law interaction to many identi- 
cal oscillators. From a purely theoretical point of view, 
the long range coupling with an appropriate normaliza- 
tion factor also allows one to interpolate between the 
extreme cases of all-to-all and nearest neighbor coupling. 
Finally, the model with power law interactions allows us 
to assess the accuracy of various analytic approximation 
schemes by comparison with large numerical simulations. 

In this paper we present a more systematic investiga- 
tion of the power law model, using analytic perturbation 
techniques analogous to spin- wave theory in magnets and 
cluster arguments following ref. [7], as well as numerical 
simulations on larger systems than in previous work. We 
study the range of interaction power laws a > 1, and the 
question of whether macroscopic synchronization may ex- 
ist for a large number of oscillators, and the nature of the 
synchronized state, as a function of a. We emphasize 
that for this range of interaction power laws, the differ- 
ent choice of normalization of the coupling strength used 
in refs. [10, 11] results only in a finite multiplicative fac- 
tor of the coupling strength, and does not change any of 
the qualitative results. Only for a < 1 is the choice of 
normalization critical to the questions being addressed. 

The outline of the paper is as follows. In Section II 
we introduce the model and the diagnostics we use to 
quantify its behavior. We then describe three approaches 
to understand the behavior. In Section III we perform 
an expansion about the aligncd-phase state analogous to 
the spin-wave approximation in magnetic systems. In 
Section IV we coarse grain the system by summing over 
blocks of oscillators, following the method used by Stro- 
gatz and Mirollo [7] in their discussion of the nearest 
neighbor model. Numerical simulations on large systems 
of up to 16384 oscillators are described in Section V. In 
Section VI we bring together the results, compare with 
previous work on the long range model [10, 11], and con- 
clude. 



II. MODEL AND DIAGNOSTICS 

In the simplest model of a population of mutually in- 
teracting oscillators with different frequencies, each oscil- 
lator is reduced to a single phase degree of freedom which 
evolves at a rate determined by its intrinsic frequency and 
its interactions with the other oscillators proportional to 
the sine of the phase differences [4] 



Yl Kij sin( 



(1) 



The all-to-all model (K^ = K/N) and the short range 
model, where only the nearest neighbors interact with 



each other (K^ = K for nearest neighbors, zero other- 
wise) have been studied in great detail [5]. 

In this paper, we consider oscillators with a power law 
coupling that varies as Kij = K/rfj , where r,j is the dis- 
tance between the oscillators at site i and j. The model 
is defined by the equations of motion 

at-1 1 

9j = uj, + K — [MOj+s ~ d i) + sin(0,- s - 6j)] ■ 

(2) 

Here, 6j (1 < j < TV) is the phase of the j th oscilla- 
tor and u>j are the corresponding intrinsic frequencies, 
assumed to be independent random variables with dis- 
tribution g(co). Without loss of generality, g(uj) can be 
chosen such that (uij) = and, for bounded distributions, 
(ujj) = 1. The parameter K sets the strength of the cou- 
pling, with a the exponent for the power-law decay of 
the interactions. We use periodic boundary conditions 
so that j + TV = j. 

We are interested in the synchronization of the oscil- 
lators to one another. A number of different criteria for 
synchronization can be introduced. 

We will use the term entrainmcnt to denote oscillators 
that are evolving with the same frequency. More pre- 
cisely, we will define oscillators i,j as entrained if there 
are no 2ir phase slips over arbitrarily long time evolution 
after initial transients have died out 
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Entrainment 

with Adij = Oi — 9j. Note that this is stricter than 
simply requiring the long-time mean frequencies u>i = 
(0i(to + T) — Oi{to))/T to be equal, since for example 
phase diffusion \A9ij(t)\ ~ t x l 2 would be consistent with 
the latter condition but not the former. A measure of 
the presence of entrainment over the whole system or 
frequency order is the Edwards- Anderson order parame- 
ter 
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For a fully entrained state all of the oscillators evolve 
with the same frequency, which will be the mean of the 
frequency distribution, and ^*ea = 1 • This is a par- 
ticularly simple state: a periodic solution (limit cycle) 
in general and a time independent solution (fixed point) 
after setting the mean frequency to zero. 

Another measure of synchronization is phase order giv- 
ing the average alignment of the phases of the oscillators. 
The degree of phase order over the system is quantified 
by the magnitude of the phase order parameter j^p/^t)! 
at any time with 
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The signal measured in an experiment that sums the 
signals from the individual oscillators is proportional to 
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Rctyphit). A state with perfectly aligned phases has 
h^p/il = 1- I 11 principle, any time dependence of H> p h 
is possible, but for the phase model we expect the syn- 
chronized motion for large N to be close to periodic. A 
fully entrained state (all oscillators entrained), for exam- 
ple, will give a periodic ^ph, but typically with \^ p h\ < 1 
since the phases will not be fully aligned. We will usu- 
ally investigate the time average magnitude (\^ P h(t)\)t 
after some time to to allow transients to decay. We also 
introduce the phase correlation function given by 

C%3 = <e^-^>, (6) 

with the average extending over time and over the lattice 
i, j for fixed separation i — j. 
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III. SPIN- WAVE ANALYSIS 

We first carry out a spin-wave (SW) type analysis of 
this model. Such an approach has been applied to the 
short range model earlier [14]. This approach studies 
the small deviations from a state of aligned phases, and 
investigates the consistency of this assumption. 



A. Preliminaries 

For large enough coupling, we might anticipate a fully 
entrained state where each oscillator evolves with the 
same frequency, and where the phase differences A8ij = 
6i — 9j between the interacting oscillators are small. The 
common frequency will be the mean of the frequency dis- 
tribution, which we have set to zero, and so the fully en- 
trained state is time independent. If the phase differences 
A9ij are small, the sine functions appearing in the inter- 
action term may be linearized. Introducing the Fourier 
modes 



(7) 



with q = 2m: /N with n integral and the sum running 
over the first Brillouin zone — ir < q ^ n, yields 



6 g =u q - K(q)O q , 
with the interaction kernel 

AT-l 



K(q) = 2K ]T 



1 — cos qs 



(8) 



(9) 
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Since K(q) is positive, each mode relaxes exponentially 
to a steady state determined by the Fourier transform of 
the random frequencies Cj q : we therefore investigate the 
properties of this steady state, which corresponds to the 
fully entrained state. Solving for the steady state of Eq. 



FIG. 1: The interaction kernel K(q) (solid) and the con- 
tinuum approximation Eq. (11) (dashed) for a = 3/2 and 
N — 4096. The inset shows the same curves on a log-log 
scale. 



(8) and using (uJ 2 ) = 1, we obtain for the mean square 



phase difference 



<I A ^| 2 > = Jf E l^)r 2 (l - cos[g(i - j)}). 



(10) 



Important issues are the behavior of (|A#y | 2 ) for large 
separations i — j, which depends on the small q terms 
in the sum in Eq. (10), and the possible divergence of 
(\A6ij\ 2 ) for any i — j due to the vanishing of K{q) for 
small q. For small q and large N, with q ^> TV -1 , we can 
evaluate K(q) in the continuum limit by replacing the 
sums by integrals with the appropriate density of states 



K{q) ~ 2K ( 
Jo 



1 — cosqs 



(11) 



The integral can be evaluated for 1 < a < 3 to give 



where 



K(q) = Kcq° 



2 sin(7ra/2)ar(— a), 



(12) 



(13) 



with r the Euler Gamma function, except at precisely 
a = 2 where K(q) = ixKq. (This is also the limit of 
the general expression for a. — > 2.) In Fig. 1, we show a 
comparison between the exact K(q) evaluated from Eq. 
(9) and the continuum approximation Eq. (12) for a = 
3/2 and N = 4096, demonstrating the accuracy of the 
approximation for A^ -1 <g<l, but over a range that 
is restricted by the finite size effects. The integral Eq. 
(11) and the sum Eq. (9) for TV — > oo both diverge for 
a < 1. 

For TV — > oo, the mean square phase difference Eq. (10) 
can be evaluated replacing the sum by an integral 
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dq\K(q)\- 2 (l-cos[q(i-j)]). (14) 
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For small q the kernel behaves as K(q) ~ g"" 1 and the 
integral Eq. (14) diverges from the small q behavior for 
a > 5/2. We interpret the divergence of (\A9ij\ 2 ) to 
signal the onset of phase slips, and the breakdown of 
the assumption of a time independent, fully entrained, 
solution. This argument therefore suggests that there is 
no fully entrained state as N —> oo for a > 5/2. 

We can also evaluate the large distance behavior of 
(\A9ij\ 2 ). For i — j — > oo the function cos[g(i — j)} is 
rapidly oscillating as a function of q and this term in Eq. 
(14) averages to zero. The remaining integral diverges, 
again because of the small q behavior, for a > 3/2, but 
is finite for a < 3/2. This implies that for a < 3/2 and 
large enough coupling strength K, the phase difference 
9i — 9j is small even for large i — j, suggesting a state 
with long range phase order, as well as entrainmcnt. On 
the other hand for 3/2<a<5/2 the argument suggests 
an entrained state without long range phase order. 

In the following sections we study the properties of 
the states in more detail, evaluating the phase correla- 
tion function and the order parameter in the spin-wave 
approximation as a function of a. 



B. Correlation Function 

The phase correlation function CV, = (cos(0j — 9j)) is 
given by 



(e 



-<|A^-|) 



(15) 



with {\A9 2 j\) obtained from Eq. (14). We now evaluate 
this expression for different ranges of a. The results are 
summarized in Eq. (26) below. 
For a < 3/2 we write Eq. (14) as 

(\A8 l3 \ 2 ) = <| AM 2 ) - - f % dq\K(q)\- 2 cos[q(i-j)}, 

71" Jo 

(16) 



with 
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dq\K(q)\ 
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the finite i — j — > oo asymptotic value. In the remain- 
ing integral in Eq. (16), the contribution from the range 
q > (i — j)^ 1 is small since the cosine term is rapidly os- 
cillating here, so that for large i—j the small q expression 
(12) can be used for K (q) and the upper limit replaced 
by oo. This gives for large r = i — j 



C r — Coo exp 



K 2 £(a)i 



3 — 2a 



1 - 



K 2 £(a)r 3 - 2a 



where 



£(a) = 7rc 2 / sin(7ra)r(3 — 2a), 



(18) 
(19) 

(20) 



with the constant c defined in Eq. (13), and 

Coo = e -<|A^ (21) 

Note that £(a) is negative for a < 3/2, so that Eq. (19) 
represents correlations growing as a power law above the 
large distance value as r is decreased. 

For a > 3/2, the full integral in Eq. (14) is dominated 
by the small q region, so that again the expression Eq. 
(12) can be used for K(q) and the upper limit replaced 
by oo. This gives 



<|A%| 2 > 



J 



2q-3 



yl ' K 2 ^(a) 
and the result for the correlation function 



(22) 



(23) 



Two special cases of note are a = 3/2 where there are 
power law correlations 



C r oc r 



(24) 



(we have not evaluated the proportionality constant), 
and a — 2 where the correlations are simply exponen- 
tial 



C r 



-r/2TT 2 K 2 



(25) 



In summary, we find the following result for the cor- 
relation function C(r) at large separations r and in the 
limit N — » oo 



C(r) 



1. 



; cxp[-l/A" 2 £(a)r 3 - 2c 



a < 1, 
1 < a < 3/2, 
a = 3/2, (26) 



= C L . 

= exp [-r 2a - 3 /K 2 £(a)] , 3/2 < 2) < 5/2, 
= exp[-r/27r 2 /v 2 ] a = 2. 



Thus as a is varied, the spin-wave approximation pre- 
dicts long range phase order for a < 3/2, and then as 
a increases the correlation function crosses over from a 
power law decay to stretched exponential, exponential, 
and then super-exponential. Similar results have been 
obtained using the spin- wave approximation for the clas- 
sical XY model [15] as the dimension is varied continu- 
ously between 1 < d < 2. In that work they were able 
to show by other methods that the stretched exponential 
prediction was an artifact of the spin-wave approxima- 
tion, and that correlations were bounded by a simple 
exponential fall off. We do not know if a similar result 
might apply to the present model, although we find some 
confirmation of the stretched exponential behavior in the 
numerical simulations presented in Section V below. 



C. Order Parameter 

An important measure of the coherence of a set of os- 
cillators is the phase order parameter Eq. (5). For an 
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FIG. 2: Phase order parameter \^ P h\ for an infinite one di- 
mensional chain calculated using the spin-wave approxima- 
tion as a function of the power-law for the decay of interac- 
tions a: solid - full calculation from Eqs. (9) and (14); dashed 
- approximation given by Eq. (28). The coupling strength is 
K = 1.0 



It is also of interest to get a more approximate es- 
timate of the order parameter in a finite system for a 
near 3/2. If the correlations remain sizable over the 
whole system we may estimate the order parameter as 
\%h(N)\ ~ y/C{N/2). (This estimate fails for large 
enough N , or when the order parameter gets small, so 
that the correlations become small over much of the 
system.) Then we approximate C(N/2) from Eq. (14): 
we use the power law form Eq. (12) for K{q) which is 
good for the small q region that dominates the inte- 
gral for a ~ 3/2; over most of the range integration 
(3ir/N < q < 7r (with /3 an 0(1) number) we argue 
that the cosine term oscillates rapidly to zero; and over 
the remainder of the range < q < (3tt/N we have 
1 — cos(qN/2) oc q 2 so that the integrand becomes small. 
This gives the estimate 



(|A^ /2 | 2 ) 



-3-2q 



ttc 2 K 2 (3 - 2a) 
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3-2a 



, (31) 



and then 



infinite system, the order parameter may be obtained 
from the asymptotic value of the correlation function 



(27) 
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2 )/2 



(32) 



Note that the expression behaves smoothly through a 
3/2, and at this value reduces to 



with Coo given by Eq. (21) and then Eq. (17). The full 
expression must be evaluated numerically. However we 
can get an analytic approximation using the approximate 
expression (12) for K(q). This should be a good approxi- 
mation for a near 3/2 where the small q range dominates 
the integral. This gives the estimate 



I* 



ph | 



Exp 



1 



T 3-2a 



27ri^ 2 [2Sin(7ra/2)ar(-a)] 2 3 - 2a 



(28) 

This approximation is compared with the result calcu- 
lated numerically using the full expression for K{q) given 
by the sum Eq. (9) for N — > oo in Fig. 2. For a — > 3/2 
Eq. (28) becomes 



ph | 



Exp 



16n 2 K 2 (3 - 2a) 



(29) 



Thus, no matter how large K is, for a close enough to 
3/2 the order parameter decreases, and goes to zero at 
a = 3/2 in the spin wave approximation. 

The magnitude of the order parameter in a finite sys- 
tem can be calculated from the correlation function 
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N 2 " £ C - 



(30) 



We have evaluated this expression, performing the dis- 
crete sums for the spin wave approximation. These re- 
sults will be compared with simulations on the dynamical 
equations in Section V below. 



(33) 



For K = 1 these expression show a very slow scaling to 
the N —> oo limits for a near 3/2: for example to achieve 
l* P h(^V)U=3/2 < 0-5 requires TV > 3 x 10 47 (setting /3 = 
I)- 



D. Self Consistency 

We might worry, if the order parameter is becoming 
small as a — > 3/2, that the spin-wave approximation 
breaks down in this vicinity due to the failure of the 
linearization of sin(#j — 9j). 

To estimate the typical size of |0, — 6j\ for the inter- 
acting oscillators, we calculate the average of the correla- 
tion function weighted by the strength of the interaction 
at that separation 



C 



ST Or 



r~ a dr 



(34) 



Evaluating the integral for the power law correlations at 
a = 3/2 given by Eq. (26), and approximating for large 
K gives 



C~ 1 



1 



4ir 2 K 2 ' 



(35) 



Thus for large enough K, the average phase correlation 
important in the interactions is close to unity, as assumed 
in the spin wave approximation, even though the order 



6 



parameter is zero. For example, C > 0.75 requires K > 

7T- 1 . 

For smaller values of a, or in a finite system, the cor- 
relations are enhanced, and so the approximation should 
be even better in these cases. 



E. General Dimension 

The results of the spin- wave calculation can be general- 
ized to d-dimensional lattices. The corresponding results 
for the correlation function are: perfect phase ordering 
with C(r) = 1 for a < d; cntrainment with long range 
phase order for a < 3d/2; a power-law fall off of the 
phase correlation function at a = 3d/2; and a cross over 
to exponential decay of correlations at a = (3d+ 1)/2 via 
stretched exponentials. 



IV. BLOCK SUMS 
A. Preliminaries 

In this section, we will analyze the long range problem 
by coarse graining the one dimensional chain into block 
oscillators. It is useful to coarse grain the chain into 
blocks by summing the equations of motion Eq. (2) for all 
6i in the block, since in this way, the internal interactions 
within a given block cancel with each other. This means 
one can look at the interaction of the block with the rest 
of the chain. In the short range model, this turns out 
to be especially useful, since the interaction of a block 
with the rest of the chain includes the surface terms only. 
In the long range model the situation is more complex, 
because the oscillators within a block interact with all 
the oscillators in the rest of the chain. In this section, for 
simplicity we use open boundary conditions rather than 
periodic ones. 

Before we discuss the results for the long range model, 
let us recapitulate the results obtained by Strogatz and 
Mirollo [7] for the nearest neighbor model. A key result 
which we will try and generalize to the long range model 
is the following. They find that it is impossible to have 
a macroscopic synchronized cluster in one dimension for 
finite values of the coupling constant. In higher dimen- 
sions, any macroscopic cluster takes the form of a sponge, 
i.e. the cluster is riddled with holes, which correspond to 
unsynchronized oscillators. 

We write the basic equations (1) in the form 



i¥=3 



sm 



(36) 



The average frequency uij of the jth oscillator is defined 
as Cjj = \im t ->oo (9 j(t) — 6j(0))/t. We define a block S as 
a contiguous segment of the chain of oscillators. For a 
synchronized block the oscillators must all have the same 
average frequency, uij = uj for all j £ S. We now sum the 



equations over a synchronized block S k of M oscillators 
entrained at frequency uj k . 

In analyzing the equation in terms of block sums we 
need the properties of two quantities: the frequency sums 
also known as the accumulated randomness; and the in- 
teraction sums. 

Summing the time averaged equations of motion (36) 
over the block S k will give on the left hand side the quan- 
tity Mwjt - Yfc(M) with 



Y k (M) 



'.r 



the accumulated randomness. We will also use 

y k {M) = M~ l Y k {M) = M" 1 £ Uj . 

jes k 



(37) 



(38) 



A key point in the arguments below is that the support 
of possible values of y k is the same as the support of 
the individual frequencies uij , but for large M the typical 
value of y k will scale as M -1 / 2 for frequency distribu- 
tions with a finite variance. For example, for a bounded 
frequency distribution, there is some probability (very 
small for large M) that each uij in the block will have 
the maximum possible value ui max , and then y k will take 
on the value ui max . On the other hand, since for large 
M the central limit theorem means that y k is given by 
a Gaussian distribution with standard deviation scaling 
as A/" 1 / 2 , the typical values of y k (those with nonzero 
probability for M — > oo) are of order M -1 / 2 . Obviously, 
similar remarks apply to Y k after including an additional 
factor of M. 

The second quantity of interest is the summed inter- 
action of the block S k with other oscillators in the chain. 
Consider first the summed interaction of this block, with 
a second block S p of size M. The largest possible inter- 
action sum, when all the phases within each block are 
aligned is 



h P = Kl 'J- 



(39) 



Using the bound 



1 rq+l/2 

£/(n)< / f(x)dx (40) 



n—p 



p-l/2 



for any function f(x) with positive curvature, we can 
bound the interaction sum for 1 < a < 2 for the two 
blocks separated by D oscillators 



I kp < K[{M + D) 2 - a + (M + Df- L 
-{M + M + Df- a - D 2 -"} 



with 



K 



K 



(2-a)(a-lY 



(41) 



(42) 
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Similarly, using the bound 

> [ 9 f(x)dx (43) 

for any monotonically decreasing function gives for 1 < 
a < 2 

J fcp >tf[(M + D) 2 ~ a + (Al + Df- a 
— (M + M + -D- 1) 2 ~ Q - (Z> + I) 2-0 ]. 

We will use various limits of these expressions. For 
example, for a block of size M interior to an infinite chain 
and interacting with the remainder of the chain we use 
D = 0, M — > oo for the interaction with the infinite chain 
in either direction, to find for 1 < a < 2 

I kp ~ 2KM 2 ~ a (45) 

(both bounds lead to the same expression in this limit). 
The essence of the block-sum arguments is to compare 
the interaction sum, scaling as M 2 ~ a , with either the 
range of possible values of the frequency sum, scaling as 
M, or the typical value, scaling as M 1 ! 2 . This will yield 
important changes of behavior at a = 1 and a = 3/2. 
For a > 2 the interaction sums between large blocks are 
independent of the sizes of the blocks, and the results es- 
sentially reduce to those of the model with nearest neigh- 
bor interactions. 



B. Impossibility of Macroscopic Clusters for Finite 
Coupling Strength and a > 1 



For a > 2 the interaction sums converge for large m, TV 
so that F a {m) = f/m with / an 0(1) constant. In the 
latter case the expression reduces to the one for the short 
range model. For a < 1, F a (m) diverges as TV — > oo. 

Now we argue that for a > 1 we can choose to suffi- 
ciently large but finite such that the probability p of Eq. 
(46) being satisfied, for block k and a given a), is less than 
unity. This follows, because there is some nonzero proba- 
bility of finding any value of yk over the support of the im- 
probability distribution, whereas the right hand side may 
be made as small as we choose by choosing to sufficiently 
large. Lemma 3.1 given in [7] makes this argument pre- 
cise. It then follows that the probability that the result 
is satisfied for all R sub-blocks is p R , and scales to zero 
as TV — > oo as 0(e~ cN ) with c some positive constant 
(remember R is 0(N)). Since the block S can be located 
at a number of different locations which is certainly less 
than TV, this means that the probability for macroscopic 
blocks satisfies P(N, K, f) < 0{Ne' cN ) which tends to 
zero as TV — > oo. 

Thus it is impossible to have a contiguous macro- 
scopic block, containing an 0(1) fraction of the oscil- 
lators locked to a common frequency, when a > 1. In 
any macroscopic segment of the chain there will always 
be (i.e. probability one as TV — » oo) some finite blocks 
of "runaway" oscillators that are desynchronized from 
their neighbors. This is the same result as in the nearest 
neighbor model: for the question of the formation of fi- 
nite blocks of unsynchronized oscillators, the power law 
interactions do not change the conclusions unless the in- 
teraction with a single oscillator can be infinite (as is the 
case for a < 1). 



In this section, we will obtain the regime of a where 
it is impossible to have a macroscopic sized synchronized 
block. More precisely, define the probability P(N, K, /) 
that there exists one or more contiguous blocks S con- 
taining at least /TV oscillators with / some finite fraction. 
Then P(N, K, f) is zero for TV ^ oo and K finite. The 
approach closely follows the one Strogatz and Mirollo [7] 
used for the short range model. 

Let us suppose that such a block S is made up of 
synchronized oscillators at a frequency a). We now di- 
vide the block into R nonoverlapping segments, Sk of 
length to each. Thus R = fN/m and k varies from 1 
to R. Note that m is an integer, sufficiently large but 
finite as TV — > oo and R is 0(N). Summing the time- 
averaged equations of motion (36) over the sub-block Sk 
and bounding the interaction sum as described above, for 
the sub-block Sk to be part of the synchronized block, we 
must have, for TV >> to and 1 < a < 2, 

\Q-y k \ <KF a (m), (46) 

with, using Eq. (41) for M = m, D = 0, M — > oo, 

2 1 

Fjm) = -— r rr. (47) 

y ' (a - 1)(2 - a) to"- 1 v ' 



C. Synchronization of Large Separated Blocks for 

a < 3/2 

For the nearest neighbor model, the result analogous to 
the one of the previous section is sufficient to show that 
for a one dimensional chain there will not be a macro- 
scopic number of synchronized oscillators for finite cou- 
pling, since the unsynchronized blocks effectively cut the 
chain into noninteracting pieces. However, for the long 
range model it is perhaps possible, even for a finite K 
and a > 1, to have a partially entrained state with a 
macroscopic number of oscillators having the same fre- 
quency. This would correspond to the system breaking 
up into disconnected blocks which however synchronize 
via the long range interaction across the unsynchronized 
oscillators. In this section we analyze the mutual interac- 
tion between two distant blocks, which are separated by 
blocks of unsynchronized oscillators, and investigate their 
possible synchronization. It is more difficult to prove the 
existence of synchronization, rather than its absence, and 
the argument we present is less rigorous than in the pre- 
vious section. 

Let us consider the two large blocks, Sk and S p , of size 
M each, where M >> 1, and a in the range 1 < a < 
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2. From the general expression Eq. (44) we obtain the 
following results. For a separation D <C M (e.g. D finite 
and M = O(N) for TV -» oo) 

Jfcp > 2if (1 - 2 1 " Q )M 2 -" (48) 

independent of the separation D with 0{D /M) correc- 
tions. Also, for blocks fc,p separated by p — fc — 1 blocks 
of size M 

I kp > c pk KM 2 - a (49) 

with c p k an O(l) number 

c pfc = 2{p-kf~ a -{p-k+l) 2 - a -{p-k-l) 2 - a . (50) 

In both cases, the lower bound on the maximum inter- 
action sum scales as M 2 ~ a . (It can be shown the upper 
bound scales in the same way). 

On the other hand the frequency sums Y k , Y m Eq. (37) 
are described, for large M, by independent Gaussian dis- 
tributions with standard deviation scaling as M 1 / 2 . The 
typical difference between the frequency sums will also 
scale as M 1 / 2 . This means that two large, fully, aligned 
blocks, each of M oscillators, separated by finite blocks 
of unsynchronized oscillators, or even by O(M) oscilla- 
tors, will typically synchronize for a < 3/2 for coupling 
strengths K > 0(M a ~ 3 / 2 ). We would expect this result 
to extend to blocks that are not fully aligned, provided 
each block has a nonzero value of the phase order pa- 
rameter , ffc, , I'p, since we would expect the interaction 
sum to be reduced by a factor of about l^fc^pl. Thus, 
for a < 3/2 the small blocks of unsynchronized oscilla- 
tors, necessarily present by the arguments of the previous 
section, do not necessarily act to break the chain into fi- 
nite lengths of synchronized oscillators, and macroscopic 
synchronization is possible. 

For a > 3/2 the typical difference of the frequency 
sums exceeds the maximum interaction sum, and so syn- 
chronization of separated blocks would not be expected 
for finite K . Comparing the scaling of the interaction 
and frequency sums with M suggests that an interac- 
tion strength scaling with system size N as N a ~ z / 2 is 
required for macroscopic synchronization for a > 3/2. 
This reduces to the result K ~ TV 1 / 2 for a = 2 where the 
block sums reduce to the nearest-neighbor results, con- 
sistent with the rigorous result of Strogatz and Mirollo 
[7] for the nearest neighbor model. 

V. NUMERICAL SIMULATIONS 

The numerical evolution of Eq. (2), has been carried 
out for five different system sizes, namely N = 256, 1024, 
4096, 8192 and 16384. Periodic boundary conditions were 
imposed. 

In order to integrate the equations of motion, we have 
used the fourth order Rungc-Kutta method with a time 
step typically At = 0.05. For the long range model with 



a system size -/V, there are C2 interaction terms. There- 
fore, evaluation of the interaction kernel using the direct 
method in real space would be an 0(N 2 ) operation at 
each time step. Instead we express the interaction as a 
convolution 

J2 K l3 sin(0< - 9-) = Im[e-^ £ K(\i - j\)e i8 *] (51) 

which can be efficiently evaluated in 0(N \ogN) opera- 
tions using fast Fourier transforms. 

The initial phases are randomly distributed between 
— 7r and 7r. The intrinsic frequencies are Gaussian ran- 
dom numbers with zero mean and unit variance. The 
first to = 200 time units are discarded as transients 
while integrating the equations. The data is then ac- 
cumulated for tf — to additional time units. We use 
tf — to = 600 for most of the results, but increase this 
to 2400 to analyze the synchronized clusters in Section 
VC. The phase order parameter has been time aver- 
aged over the entire range tf — to, while the Edwards- 
Anderson order parameter has been calculated as f ^ = 
I X^li eW&ft-t'VoVl/N. The data for the order param- 
eters and the correlation functions have also been aver- 
aged over 30 different initial configurations, by repeating 
with different seeds for the random number generator. 
We present the results for a coupling strength K = 1. 

A. Order Parameters 




a 



FIG. 3: The phase order parameter \^ P h\ as a function of 
the power law a of the interaction for five different system 
sizes. The strength of the coupling is K = 1. The system 
size increases from the upper curve to the lower: squares - 
256; circles - 1024; triangles - 4096; pluses - 8192; diamonds - 
16384. 

We show the magnitude of the phase order parameter 
Eq. (5) as a function of the power law a of the interac- 
tion, for coupling strength K = 1 and for the five different 
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a 

FIG. 4: Estimate of the phase order parameter \9 p h\ (solid 
curves) as a function of the power law a of the interaction 
based on Eqs. (31) and (32) for the same five system sizes 
used in Fig. 3, increasing from the upper curve to the lower. 
The strength of the coupling is K = 1 and a value of j3 = 0.5 
was used. The spin wave prediction for an infinite system size 
is shown for reference (dashed curve). 



system sizes in Figs. 3. The order parameter decreases 
rapidly for a value of a that decreases with increasing 
system size, with the half-height \^ p h\ = 0.5 value occur- 
ring at about a = 1.6 for the largest system size used. As 
we have seen in Section III C we expect strong finite size 
effects for a near 3/2, so that it is hard to extrapolate the 
numerical data to the infinite size limit. For comparison 
we show in Fig. 4 the finite size estimate of the order 
parameter from Eqs. (32) and (31). We do not expect 
a quantitative agreement for any a due to the crude ap- 
proximations made, and the behavior for small \^ P h\ is 
not correct as explained in Section IIIC, but the overall 
trends are quite similar. 

A comparison between the full spin- wave predictions, 
performing the discrete sums without approximations, 
and the results from the simulations of the time evolu- 
tion, both for system size 8192, is shown in Fig. 5. The 
plot shows quantitative agreement between the simula- 
tions and the spin wave sums for a < 1.4. As a ap- 
proaches 1.5 and the order parameter decreases, there is 
increasing disagreement, as would be expected since the 
theory is an expansion assuming well aligned phases. 

Although the strong finite size effects preclude reliably 
extrapolating to infinite system sizes, the numerical simu- 
lation results appear to be consistent with the prediction 
of the spin wave theory, that phase ordering and entrain- 
ment is possible for a < 3/2, but that phase order is not 
possible for N — > oo for a > 3/2. 

The Edwards- Anderson order parameter evaluated 
from the simulations for the same parameter values is 
shown in Fig. 6. For all system sizes used and coupling 
strength K ~ 1, ^ ea remains unity up to a = 3/2, 
showing perfect entrainment, suggesting a single block 




0.0 -h . , . 1 . i JJJ -' 7 
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a 



FIG. 5: Comparison between the spin-wave sum (dashed line) 
and the simulation results (squares joined by solid line) for the 
phase order parameter \^ p h\ as a function of the power law a 
of the interactions in a system of size 8192 and for coupling 
strength K = 1. The inset shows the same comparison over 
a smaller range 




a 



FIG. 6: The Edwards- Anderson order parameter vPea as a 
function of the power law a of the interaction for five different 
system sizes as in Fig. 3 and coupling strength K = 1. 



of all the oscillators evolving at the mean frequency of 
the distribution for the system sizes used. For a > 3/2 
the order parameter decreases, showing that the system 
has broken up into more than one synchronized cluster. 
Qualitatively, the fall off above a = 3/2 is similar to the 
fall off of the phase order parameter. Thus we do not find 
stronger correlations in the oscillator frequencies than in 
the phases. This does not appear to be consistent with 
the predictions from the spin- wave theory of a state with 
entrained oscillators but without long range phase order 
for the range 3/2 < a < 5/2. On the other hand the 
result is consistent with the block sum arguments which 
show that for a < 3/2 blocks of oscillators may synchro- 
nize across regions of unsynchronized oscillators, whereas 
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this does not typically happen for a > 3/2. 



B. Correlation Function 
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FIG. 7: Comparison between the spin-wave sum and the sim- 
ulation results for the correlation function C r as a function 
of separation r for six values of the interaction power law ex- 
ponent a and a system size of 4096 and coupling strength 
K = 1: dark lines - spin wave prediction; light points - simu- 
lation results. 

As a further test of the results of the spin wave ap- 
proach we show in Fig. 7 a comparison between the 
phase correlation function Eq. (15) evaluated from the 
spin wave discrete sums and the full numerical simula- 
tions for six values of a in the range 1 < a < 3/2 for a 
system size 4096. The plots show quantitative agreement 
for a < 1.3. As a approaches 1.5, the difference grows, 
which is consistent with what we saw earlier in Fig. 6. 

For 3/2 < a < 2 the spin wave analysis predicts a 
phase correlation function with a stretched exponential 
decay, see Eq. (26). Based on these predictions, we fit 
the correlation function obtained from the simulations in 
this range of a to a function 



C(r) 



— br c 



(52) 



for a — 1.6, a = 1.7, and a = 1.8 and system size 8192, 
see Fig. 8. An exponential correlation would be a straight 
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FIG. 8: Phase correlation function as a function of separation 
for system size 8192: light points - simulations; dark lines - 
stretched exponential fits, Eq. (52). The values of the inter- 
action power law a are, top to bottom, a = 1.6, a = 1.7, and 
a = 1.8. 





b(Numerical) 


b(SW) 


c(Numerical) 


c(SW) 


a = 1.6 
a = 1.7 
a = 1.8 


0.094 ± 0.001 
0.151 ±0.001 
0.151 ±0.001 


0.093 
0.064 
0.055 


0.330 ± 0.001 
0.456 ± 0.001 
0.583 ± 0.002 


0.2 
0.4 
0.6 



TABLE I: Parameters of the stretched exponential correlation 
function Eq. (52) from fits to numerical results for system size 
8192 and three values of a > 1.5 compared with predictions 
from spin-wave theory for the infinite system. 



line on these log-linear plots, and clearly does not fit the 
data. The fit parameters for the stretched exponential 
and the predictions of the spin wave theory are shown in 
Table I. The agreement between the exponents is rea- 
sonably good for a = 1.7 and a = 1.8. The deviation 
from the predictions for a = 1.6 can be ascribed to finite 
size effects, since the correlations tend to a finite value for 
large r corresponding to separations of N/2 for this value 
of a, and so should not be compared with the theoretical 
predictions for an infinite size system. In the next sec- 
tion we discuss the size of the synchronized clusters. The 
range over which the stretched exponential fit is good in 
Fig. 8 is comparable to the maximum cluster size. 



C. Clusters 

In this section we present results from our simulations 
for the size of synchronized clusters as a function of the 
power law of the interaction. For long range interactions 
both contiguous blocks, and disjoint blocks entrained 
through the long range interaction across unentraincd os- 
cillators, are of interest. 

We identify an entrained cluster from the simulations 
as a set of oscillators that are phase locked: over the 
time of the simulation no oscillator phase undergoes slips 
(changes of about ±27r) with respect to the mean phase. 
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In a simulation over a time T this is equivalent to the fre- 
quency being within 2ir/T of the mean frequency of the 
cluster. The expected frequency difference between large 
but distinct clusters of size about L is of order L -1 / 2 . 
Thus the simulation time should exceed 2-kL 1 / 2 . We com- 
pute the phase-winding number, n w , for every oscillator 
along the chain. The phase- winding number is calculated 
as 



lim* 



,(&(*) -0i(t o )) 



4tt 



(53) 



where [x] denotes the nearest integer to x. 
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FIG. 9: Winding numbers over a portion of the system for 
a = 1.7 and K = 1 and a system size of 8192. 

An example of the raw data of winding numbers, com- 
puted over a run time T = t — 1 of 2400, is shown in Fig. 
9 for a system size N = 8192, coupling strength K = 1 
and interaction power law a = 1.7, Only a portion of the 
full system is shown. Note in particular blocks with the 
same winding numbers entrained across oscillators with 
different winding numbers, a novel consequence of the 
long range interaction. 

Histograms of the phase winding number for the same 
system are shown in Fig. 10 for four values of a. We use 
the total number of oscillators with the same winding 
number, shown by the bar length, to define the overall 
cluster size, and this is divided up into the individual 
contiguous blocks (containing no oscillators with differ- 
ent winding numbers) given by the lengths between the 
points on the bars. For clarity, only a restricted range 
of winding number is shown: there are additional small 
clusters with more distant winding numbers outside the 
range plotted. For a = 1.6 almost all the oscillators have 
the same winding number n w = 0, so that the cluster of 
entrained oscillators spans the whole system, with only a 
few oscillator of different winding numbers breaking the 
global cluster into four smaller contiguous blocks. As a 
increases, more clusters of smaller size develop. In each 
case a number of contiguous blocks join to form a large 
cluster with entrainmcnt across uncntraincd oscillators. 
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FIG. 10: Phase winding numbers Tin) ELS Qb function of cluster 
size for a system of size 8196 and interaction strength K = 1. 
The lengths of the bars define the total number of oscillators 
with the same winding number, and the lengths between the 
points on the bars give the sizes of the individual contiguous 
blocks making up the whole. The four plots are for different 
interaction power laws: (a) a = 1.6; (b) a = 1.7; (c) a = 1.8; 
(d) a = 1.9. Only a restricted range of winding number is 
shown. 



The full distribution of contiguous block sizes is shown 
in Fig. 11. This is a plot of an ordered list of the sizes 
of contiguous blocks. A log-linear plot of the same data 
shows a good fit to an exponential fall off for small block 
sizes (the large block number end of the plot). 



VI. CONCLUSIONS 

We have studied the synchronization of oscillators de- 
scribed by a phase only model with interactions falling 
off with separation r as a power law r~ a using an expan- 
sion about the aligned phase state (spin- wave method), 
arguments summing the equations of motion of blocks 
of oscillators (block-sum method) and numerical simu- 
lations on systems of up to 16384 oscillators. We have 
focusscd on the range a > 1, since previous work [11] has 
looked at a < 1 in some detail. 

For 1 < a < 3/2 we find results consistent with macro- 
scopic entrainment and long range phase order for large 
enough coupling strengths. The spin-wave type analysis, 
based on the assumption of a time independent solution 
(fully entrained state) and an expansion of the interaction 
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FIG. 11: Ordered plot of contiguous block sizes for four values 
of the power law a: squares - a = 1.6; circles - a = 1.7; 
triangles - a = 1.8; diamonds - a = 1.9. The ordinate is 
the size of a contiguous entrained block, and the abscissa the 
index number in the list of blocks ordered by size. Other 
parameters are as in Fig. 10. 



term linearly in 6i — 9j, predicts a state with long-range 
phase order and a nonzero phase order parameter. For 
large enough coupling strength K the mean square phase 
deviation ((8i — 9j) 2 ) is small, and the average phase cor- 
relations weighted by the power law interaction are close 
to unity, so that the linear expansion of the nonlinear 
interaction function is a good approximation on average. 
The block sum argument shows that despite the long 
range interaction, the coupling of a finite block to the 
rest of the chain remains bounded above by some finite 
value, and there is a nonzero probability of finding a finite 
block of oscillators with frequencies sufficiently far from 
the mean that they are not synchronized to the rest of 
the chain for finite coupling. This argument shows that 
for N —> oo, there are no macroscopic (O(N)) contiguous 
blocks of synchronized oscillators for any finite K, so that 
the assumption of a time independent solution as made 
in the spin- wave approach is not correct [16]. However, 
for a < 3/2, the interaction is sufficiently long range that 
large blocks of oscillators are likely to synchronize across 
the unsynchronized oscillators (see Fig. 9 for examples 
from the simulations), leading to the entrainmcnt of a 
finite fraction of the oscillators, long range phase correla- 
tions, and a nonzero order parameter for sufficiently large 
K, even for N — » oo. These results follow the predictions 
of the spin- wave theory, although the finite blocks of un- 
synchronized oscillators will reduce the order parameter 
and phase correlations below the value predicted by spin 
wave theory, as seen in the comparison of the simulations 
with the spin-wave predictions. 

For 3/2 < a < 5/2 the spin- wave approach predicts a 
fully entrained state, but with no long range phase order, 
although for a < 2 the phase correlations are predicted to 
be of strctched-exponential form. However the block sum 



method shows that finite unsynchronized blocks again ex- 
ist, and now large blocks of oscillators will typically not 
synchronize across the unsynchronized oscillators. Thus 
for finite coupling strength K and a > 3/2 we expect 
no macroscopic entrainment (no finite fraction of oscilla- 
tors at the same frequency for N —> oo). The results of 
the numerical simulations show the Edwards- Anderson 
order parameter which measures frequency entrainmcnt, 
decreasing as a increases above 3/2 in a way that is 
broadly similar to the phase order parameter, consistent 
with the picture that the unsynchronized blocks disrupt 
both the phase and frequency correlations. The simula- 
tions show results consistent with the spin-wave predic- 
tions of a stretched exponential decay of correlations up 
to a distance comparable with the largest cluster size. 

Rogers and Willc concluded in their paper that the 
critical interaction exponent a c such that the oscillators 
do not synchronize for a > a c even for very large coupling 
strengths is a c ~ 2. Our results suggest a lower critical 
value of a c = 3/2 and our numerics on larger systems 
than used in ref. [10] approach this value. The diagnos- 
tic used by Rogers and Wille was the average plateau 
size as a fraction of the system size. In their simula- 
tions they found this quantity to switch quite rapidly as 
a function of increasing a for reasonably large K from 
unity to close to zero. The plateaus were defined as con- 
tiguous blocks of oscillators with the same frequency, and 
so oscillators synchronized across unsynchronized blocks 
were not counted as in the same plateau. This means 
that their diagnostic does not detect long range synchro- 
nization occurring through this mechanism. However, we 
believe the main reason that their value of a c is greater 
than the value 3/2 that we propose is due to the strong 
finite size effects for a near 3/2, so that for the range of 
sizes they used, too large a value of a is needed for desyn- 
chronization to appear, and their extrapolation scheme 
to large N was not adequate. 

Marodi et al. have also looked at this problem numeri- 
cally for sizes up to 1000 in one dimension (as well as two 
dimensions). The main focus of their work was a < 1, 
where complete entrainmcnt occurs for TV — > oo for the 
scaling of the coupling constant they use. They do not 
attempt to identify a critical value of a > 1 above which 
partial synchronization is no longer possible. Their nu- 
merics on system sizes up to 1000 and for K = 1 show 
the phase order parameter decreasing for a close to 3/2 
— in fact closer than we find for these system sizes (com- 
pare their Fig. 2 with our Fig. 3) perhaps because of the 
open boundary conditions they use. They also remark 
that the order parameter approaches a steady value for 
a > 2.5, but we believe this value tends to zero for large 
N, which is consistent with the trends in their numerics. 
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